\(N = 3154\) healthy young men aged 39–59 from the San Francisco area were assessed for their personality type. All were free from coronary heart disease at the start of the research. Eight and a half years later change in this situation was recorded.
behave: behavior type which is a factor with levels A1 A2 B3 B4
cigs: number of cigarettes smoked per day
dibep: behavior type a factor with levels A (Agressive) B (Passive)
chd: coronary heat disease developed is a factor with levels no yes
typechd: type of coronary heart disease is a factor with levels angina infdeath none silent
timechd: time of CHD event or end of follow-up
arcus: arcus senilis is a factor with levels absent present
For now, we’ll focus on 3 variables
Show R code
summary(wcgs[, c("chd", "height", "cigs")])
chd height cigs
no :2897 Min. :60.00 Min. : 0.0
yes: 257 1st Qu.:68.00 1st Qu.: 0.0
Median :70.00 Median : 0.0
Mean :69.78 Mean :11.6
3rd Qu.:72.00 3rd Qu.:20.0
Max. :78.00 Max. :99.0
Anything interesting stick out?
Visualizing discrete data
Pie chart of response 🤮
Show R code
# Construct a pie chart of the (binary) response; I'm not a fan of pie charts in generalptab <-prop.table(table(wcgs$chd)) # convert frequencies to proportionsptab # inspect output
no yes
0.91851617 0.08148383
Show R code
pie(prop.table(table(wcgs$chd)),main ="Pie chart of Coronary Heart Disease")
Visualizing discrete data
Bar charts tend to be more effective
Show R code
barplot(ptab, las =1, col ="forestgreen")
Visualizing discrete data
Mosaic plot showing relationship between cigs and chd; not incredibly useful IMO unless both variables are categorical
Show R code
plot(chd ~ cigs, data = wcgs)
Visualizing discrete data
Nonparametric density plot of height by chd status using lattice graphics:
Show R code
library(lattice)densityplot(~ height, groups = chd, data = wcgs, auto.key =TRUE)
Visualizing discrete data
Boxplot of cigs vs. chd status
Show R code
plot(cigs ~ chd, data = wcgs, col =c(2, 3))
Visualizing discrete data
Boxplot of cigs vs. chd status
Show R code
plot(cigs ~ chd, data = wcgs, col =c(2, 3))
Visualizing discrete data
Boxplot of height vs. chd status with notches
Show R code
plot(height ~ chd, data = wcgs, col =c(2, 3), notch =TRUE)
Visualizing discrete data (detour)
Decision trees🌲🌴🌳 are immensely useful for exploring new data sets!
plot(partykit::ctree(chd ~ height + cigs, data = wcgs))
Observations so far…
It seems that the cigs is positively associated with the binary response chd
It is not clear how, if at all, height is associated with chd
Question: how can we build a model to examine these potential associations?
Linear models
Recall that in linear regression we model the conditional mean response as a linear function in some fixed, but known parameters \(\boldsymbol{\beta}\):
If \(Y\) is a binary random variable, then what is \(E\left(Y|\boldsymbol{x}\right)\)?
The linear probability (LP) model
It turns out that \(E\left(Y|\boldsymbol{x}\right) = P\left(Y = 1|\boldsymbol{x}\right)\). The LP model assumes that \[
P\left(Y = 1|\boldsymbol{x}\right) = \beta_0 + \beta_1x_1 + \dots \beta_px_p = \boldsymbol{\beta}^\top\boldsymbol{x}
\]
Is this reasonable?
The LP model: chd vs. cigs
Show R code
y <-ifelse(wcgs$chd =="yes", 1, 0)summary(fit <-lm(y ~ cigs, data = wcgs))
Call:
lm(formula = y ~ cigs, data = wcgs)
Residuals:
Min 1Q Median 3Q Max
-0.25268 -0.09794 -0.05876 -0.05876 0.94124
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0587608 0.0062041 9.471 < 2e-16 ***
cigs 0.0019588 0.0003339 5.867 4.91e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2722 on 3152 degrees of freedom
Multiple R-squared: 0.0108, Adjusted R-squared: 0.01049
F-statistic: 34.42 on 1 and 3152 DF, p-value: 4.91e-09
Are the standard errors here appropriate? Why/why not?
The LP model: chd vs. cigs
Show R code
plot(y ~ cigs, data = wcgs, ylim =c(0, 1), las =1)abline(fit, col =2)
The logistic regression (LR) model
Assume \(Y \sim \mathrm{Bernoulli}\left(p\right)\), where \[
p = p\left(\boldsymbol{x}\right) = P\left(Y = 1|\boldsymbol{x}\right)
\] and \[
\mathrm{logit}\left(p\right) = \log\left(\frac{p}{1-p}\right) = \boldsymbol{\beta}^\top\boldsymbol{x}
\] In other words, LR models the logit of the mean response as a linear function in \(\boldsymbol{\beta}\); we’ll refer to the term \(\eta = \boldsymbol{\beta}^\top\boldsymbol{x}\) as the linear predictor. Why does this make more sense?
The logistic regression (LR) model
Can always solve for \(p\) to get predictions on the raw probability scale (Homework 2😄): \[
p\left(\boldsymbol{x}\right) = \frac{\exp\left(\boldsymbol{\beta}^\top\boldsymbol{x}\right)}{1 + \exp\left(\boldsymbol{\beta}^\top\boldsymbol{x}\right)}
\]
Note how the LR model is nonlinear in \(p\)!
Fitting an LR model in R
Use the glm() function instead of lm()
GLM stands for generalized linear model, which includes the LR and ordinary linear regression models as special cases
Many (but not all) of the models we’ll discuss in throughout this course belong to the class of GLMs
Note how we have to specifcy the family argument! (see ?glm for details)
summary(fit.lr <-glm(chd ~ cigs, data = wcgs, family = binomial))
Call:
glm(formula = chd ~ cigs, family = binomial, data = wcgs)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.9958 -0.4418 -0.3534 -0.3534 2.3684
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.742160 0.092111 -29.770 < 2e-16 ***
cigs 0.023220 0.004042 5.744 9.22e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1781.2 on 3153 degrees of freedom
Residual deviance: 1750.0 on 3152 degrees of freedom
AIC: 1754
Number of Fisher Scoring iterations: 5
Fitting an LR model in R
Let’s compare the fitted LR and LP models:
Show R code
prob <-predict(fit.lr, newdata =data.frame(cigs =0:99), type ="response")plot(y ~ cigs, data = wcgs, ylab ="chd", las =1, xlim =c(0, 99))lines(0:99, y = prob, col =2)abline(fit, col =3, lty =2)legend("topright", legend =c("LR fit", "LP fit"), lty =c(1, 2), col =c(2, 3))
Fitting an LR model in R
Let’s compare the fitted LR and LP models:
Show R code
prob <-predict(fit.lr, newdata =data.frame(cigs =0:999), type ="response")plot(y ~ cigs, data = wcgs, ylab ="chd", las =1, xlim =c(0, 999))lines(0:999, y = prob, col =2)abline(fit, col =3, lty =2)legend("topright", legend =c("LR fit", "LP fit"), lty =c(1, 2), col =c(2, 3))
Fitting an LR model in R
Now let’s include an additional predictor (i.e., height):
Show R code
summary(fit.lr2 <-glm(chd ~ cigs + height, data = wcgs, family = binomial))
Call:
glm(formula = chd ~ cigs + height, family = binomial, data = wcgs)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.0041 -0.4425 -0.3630 -0.3499 2.4357
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.50161 1.84186 -2.444 0.0145 *
cigs 0.02313 0.00404 5.724 1.04e-08 ***
height 0.02521 0.02633 0.957 0.3383
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1781.2 on 3153 degrees of freedom
Residual deviance: 1749.0 on 3151 degrees of freedom
AIC: 1755
Number of Fisher Scoring iterations: 5
Interpreting LR coefficients
Let \(p = P\left(Y = 1\right)\) and \(1 - p = P\left(Y = 0\right)\)
The odds of \(Y = 1\) occuring is defined as \(p / \left(1 - p\right)\)
For a fair coin 🪙, the probability of getting tails is \(p = 0.5\). Therefore, the odds of getting tails vs. heads is \(p / (1 - p) = 0.5 / 0.5 = 1\). (We might also say the odds of getting tails is “1 to 1”).
For a fair die 🎲, what are the odds of rolling a 2?
Interpreting LR coefficients
The logit models the log odds of success (i.e., \(Y = 1|\boldsymbol{x}\)) \[
\log\left(\mathrm{odds}\right) = \log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \dots \beta_px_p
\] Exponentiating both, we get \[
\mathrm{odds} = \frac{p}{1-p} = \exp{\left(\beta_0\right)}\times\exp{\left(\beta_1x_1\right)}\times\exp{\left(\beta_2x_2\right)}\times\dots\times\exp{\left(\beta_px_p\right)}
\]
In the LR model, \(\beta_i\) represents the change in the log odds when \(x_i\) increases by one unit (all else held constant)
In the LR model, \(\exp\left(\beta_i\right)\) represents the multiplicative increase in the odds when \(x_i\) increases by one unit (all else held constant)
CANNOT interpret the coefficients in terms of \(p\) directly…effect plots to the rescue 🛟!!
All have there own assumptions and drawbacks; typically, similar in shape when the model is additive in nature (i.e., no interaction effects)
Effect plots
The plotmo library is an “easy button” for quick and dirty effect plots (other variables are held fixed at their median or most frequent value) and supports a wide range of models
Maximizing the log-likelihood is equivalent to minimizing the negative log-likelihood (which is more convenient):
Show R code
# Response (as a binary 0/1 variable)y <-ifelse(wcgs$chd =="yes", 1, 0)# Model matrix; includes a column for the intercept by defaultX <-model.matrix(~ cigs + height, data = wcgs)# Function to compute the negative log-liklihood (as a function of the betas)nll <-function(beta) { lp <- X %*% beta # linear predictor; same as b0 + b1*x1 + b2*x2 + ...-sum(y * lp -log(1+exp(lp)))}# Use general optimization; would be better to use gradient and hessian info (e.g., first and second derivative info)lp.fit <-lm(y ~ cigs + height, data = wcgs)optim(coef(lp.fit), fn = nll, control =list("maxit"=9999, "reltol"=1e-20))
$par
(Intercept) cigs height
-4.50160987 0.02312741 0.02520773
$value
[1] 874.5246
$counts
function gradient
679 NA
$convergence
[1] 0
$message
NULL
Extending the general F-test
In LR, we move from a general linear F-test to a likelihood ratio test based on the \(\chi^2\) distribution:
Show R code
anova(fit.lr, fit.lr2, test ="Chi")
Confidence intervals
Lot’s of ways to do this (some bad, some better). R defaults to using a profile likelihood method:
